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Abstract 

We apply the concept of effective order to strong stability preserving (SSP) explicit Runge-Kutta 
methods. Relative to classical Runge-Kutta methods, effective order methods are designed to satisfy 
a relaxed set of order conditions, but yield higher order accuracy when composed with special starting 
and stopping methods. The relaxed order conditions allow for greater freedom in the design of effective 
order methods. We show that this allows the construction of four-stage SSP methods with effective 
order four (such methods cannot have classical order four). However, we also prove that effective 
order five methods — like classical order five methods — require the use of non-positive weights and so 
cannot be SSP. By numerical optimization, we construct explicit SSP Runge-Kutta methods up to 
effective order four and establish the optimality of many of them. Numerical experiments demonstrate 
the validity of these methods in practice. 

1 Introduction 

Solutions of non-linear hyperbolic partial differential equations (PDEs) may contain discontinuities even 
when the initial conditions are smooth. The challenge for the numerical solution of these systems is 
twofold. It is desirable that the approximation be of high accuracy in regions where the solution is 
smooth and that the discontinuities be captured without exhibiting any oscillations or overshoots. 

There has been considerable effort to develop, for these and other problem classes, numerical methods 
which are strongly stable. Many of these numerical methods are based on a method-of-lines approach 
where the problem is first discretized in space to yield a system of ODEs. The spatial discretization is 
often chosen to ensure certain strong stability properties of the original PDE problem (e.g., max- norm 
monotonicity, total variation boundedness, positivity, etc.) are preserved when coupled with first-order 
forward Euler time integration. Strong stability preserving (SSP) time discretizations (formerly TVD 
discretizations [10]) are high-order time discretizations that guarantee the same stability preservation, 
with a possibly different step-size restriction. 

We examine the SSP properties of explicit Runge-Kutta methods of effective order. Effective order 
methods use special starting and stopping procedures in such a way that the method can achieve an 
order of accuracy higher than its classical design order. This allows the construction of high-order SSP 
Runge-Kutta schemes by using low-order SSP Runge-Kutta methods. 
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Explicit SSP Runge-Kutta methods of classical order four require at least five stages [10]. In contrast, 
we construct explicit SSP Runge-Kutta methods of effective order four with only four stages. Following 
this result, we had hoped to overcome the fifth-order barrier for explicit SSP Runge-Kutta methods [23]; 
instead we prove that the barrier also holds for SSP methods of effective order five. Most of the methods 
we find are optimal, as they achieve a certain theoretical upper bound on the SSP coefficients that is 
obtained by considering only linear problems [19]. 

The rest of the paper is organized as follows. Section 2 reviews Runge-Kutta methods and the concept 
of strong stability preserving methods. Section 3 presents a brief overview of the algebraic representation 
of Runge-Kutta methods, following Butcher [4]. This includes the concept of effective order and a list 
of effective order conditions. Section 4 proves an order barrier for effective order methods with strictly 
positive weights, a consequence of which is the non-existence of explicit SSP Runge-Kutta methods of 
effective order five. Section 5 presents effective order SSPRK methods found by numerical search, some of 
which are established as optimal. Starting and stopping methods are also discussed. The paper concludes 
with numerical experiments in Section 6 and conclusions in Section 7. 

2 Strong stability preserving Runge-Kutta methods 

Strong stability preserving (SSP) time-stepping methods were originally introduced for time integration 
of systems of hyperbolic conservation laws [25] 

U t + V • f(U) = 0, (2.1) 

with appropriate initial and boundary conditions. A spatial discretization gives the system of ODEs 

u'(t)=F(t,u(t)), (2.2) 

where u is a vector of continuous-in-time grid values approximating the solution U at discrete grid points. 
Of course, (2.2) can arise in many ways and F need not necessarily represent a spatial discretization. In 
any case, a time discretization then produces a sequence of solutions u n « u(t n ). In this work we study 
explicit Runge-Kutta time discretizations. An explicit s-stage Runge-Kutta method takes the form 

5 

u n+l = u n +At J2 b t F{t n + CiAt, Yi), 

i 

where 

Yi = u n + At a tj F(t n + Cj At, Y 3 ) 

3 

and Ci — X)j=i a ij ■ The accuracy and stability of the method depend on the coefficients (A, b, c) [4] . 

In some cases, the solutions of hyperbolic conservation laws satisfy a monotonicity property. For 
example, if (2.1) is scalar then solutions are monotonic in the total variation semi-norm [15]. For this 
reason, many popular spatial discretizations are designed such that, for a suitable class of problems, the 
solution u in (2.2) computed with the forward Euler scheme is non-increasing (in time) in some norm, 
semi- norm, or convex functional; i.e., 

+ AtF(t,u)\\ < ||u||, for all u and for < At < At FE . (2.3) 

Note that AipE is a property of F (and is independent of it). If this is the case, then an SSP method 
also generates a solution whose norm is non-increasing in time, under a modified time-step restriction. 

Definition 2.1 (Strong Stability Preserving). A Runge-Kutta method is said to be strong stability 
preserving with SSP coefficient C>0if, whenever the forward Euler condition (2.3) holds and 

0<At< CAt FE , 
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the Runge-Kutta method generates a monotonic sequence of solution values u n satisfying 

\\u n+1 \\ < \\u n \\. 

The SSP coefficient C is a property of the particular time-stepping method and quantifies the allowable 
time step size relative to that of the forward Euler method. Generally we want the SSP coefficient to be 
as large as possible for efficiency. To allow a fair comparison of different explicit methods, we consider 
the effective SSP coefficient 

C c ff = — • 
s 

Note that the use of the word effective here is unrelated to the concept of effective order introduced in 
Section 3. 

2.1 Optimal SSP schemes 

We say that an SSP Runge-Kutta method is optimal if it has the largest possible SSP coefficient for a 
given order and a given number of stages. The search for these optimal methods was originally based 
on expressing the Runge-Kutta method as combinations of forward Euler steps (the Shu-Osher form) 
and solving a non-linear optimization problem [10,11,22,24,26,27]. However, the SSP coefficient is 
related to the radius of absolute monotonicity [20] and, for irreducible Runge-Kutta methods, the two 
are equivalent [7, 14] . This gives a simplified algebraic characterization of the SSP coefficient [8] ; it is the 
maximum value of r such that the following conditions hold: 

K{I + rA)- l >Q (2.4a) 
e s+x -rK(I + rA)~ 1 e a > 0. (2.4b) 

Here 




while e s denotes the vector of ones of length s. The inequalities are understood component- wise. 

The optimization problem of finding optimal SSP Runge-Kutta methods can thus be written as 
follows: 

max r subject to (2.4) and $(K) = 0. (2.5) 

A,b,r 

Here Q(K) represents the order conditions. 

Following [15,18], we will numerically solve the optimization problem (2.5) to find optimal effective 
order explicit SSP Runge-Kutta methods. However, we first need to define the order conditions $(K) 
for methods of effective order. This is discussed in the next section. 

3 Effective order Runge-Kutta methods 

The effective order of a Runge-Kutta method is defined in an abstract algebraic context introduced by 
Butcher [1] and developed further in [2,3,5, 13] and others. In this section we follow [4], reviewing the 
fundamental concepts of this representation which are then used to define effective order methods and 
their order conditions. 

3.1 The Runge-Kutta group 

Runge-Kutta methods can be viewed as elements of an algebraic group in which the product corresponds 
to the composition of two methods. Let G be the group of all real- valued maps on the set of rooted trees. 
Each function a £ G corresponds to an equivalence class of Runge-Kutta methods and maps trees to 
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Table 3.1: Elementary weights of trees up to order five for a Runge— Kutta method (A,b,c). Here C 
is a diagonal matrix with components c$ = X)j=i a ij an< ^ exponents of vectors represent component 
exponentiation. 

specific algebraic expressions in the coefficients of a Runge-Kutta method, known as elementary weights. 
Two Runge-Kutta methods belong to the same equivalence class if they have the same elementary weight 
values. 

For every function a £ G we write the values of the elementary weights as cti — a(ij) for trees 
ti indexed by integer i. The vector of these values a*, i = 0, 1,2, . . . is referred to as the B-series of 
the corresponding Runge-Kutta method and is related to Taylor expansions of the numerical solution 
given by the Runge-Kutta method [4,13]. By convention a(io) = 1, where to denotes the empty tree. 
Table 3.1 lists these expressions for trees of up to degree five; a general recursive formula can be found 
in [4, Definition 312A]. 

Let a, (3 S G correspond to Runge-Kutta methods Mi and respectively. A multiplicative group 
operation a/3 can be defined by partitioning the input tree and computing over the resulting forest [4]. 
This product a/3 is related to the application of method M\ followed by method M2; we denote the 
resulting method by M2M1. 1 The product is defined by 

(a/3)(i) = W JJ a(v){3(w)\ (3.1) 

where w < t indicates a subtree of t which includes the root of t and w\t indicates the forest induced 
by removing w from t [4]. Multiplicity in choosing w must also be accounted for. The following example 
shows how to compute this product for one particular tree. 

Example 3.1. Table 3.2 shows the partition of the five-vertex tree t\\ to all possible rooted subtrees. 
Based on this partition, we apply (3.1) to find that the product of two functions in G on tree t\\ is given 
by (a/3) (in) = «n + aia 3j 8i + (a? + a 3 )/3 2 + af/3 3 + 2a^/3 4 + 2ai/3 6 + ai/3 7 + /3 n . 

3.2 Algebraic interpretation of order 

If two Runge-Kutta methods correspond to the same element in G, then they are essentially the same 
method (up to reducibility). However, the definition of equivalence of methods is overly restrictive for 
practical purposes. A weaker condition is established if we discuss equivalence of methods up to a 
particular order of accuracy. 

1 We write M^M\ to mean the application of method M\ followed by the application of method M.% (following matrix 
and operator ordering convention) but when referring to products of elements of the Runge-Kutta group we use the reverse 
ordering (a/3) to match the convention in [4]. 
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Table 3.2: Partitions of the tree t\\ to all possible subtrees w and the corresponding forests t\x \ w. 
Multiplicity is indicated with (x2). 
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Definition 3.2. Two Runge-Kutta methods Mi and M2, are equivalent up to order p if their correspond- 
ing elements in G, a and j5, satisfy 

a{t) = (3{t), for every tree t with r(t) < p, 

where r(t) denotes the order of the tree (number of vertices). We denote this equivalence relation by 
Mi ~ M 2 . 

V 

In this sense, methods have inverses: the product of a -1 and a must match the identity method up 
to order p. Classical order follows from comparing a method with the special group element E 6 G which 
advances the exact solution by one step. All of this can be made considerably more precise using quotient 
groups of G [4]. 

Example 3.3. Consider the forward Euler method u n+1 = u" + AtF(u n ). To find an inverse, we seek 
a method that undoes the work of this method, recovering u n from u n+1 ; one approach is to solve for 
u n , obtaining the backward Euler method with a time-step of —At. Alternatively, let a e G correspond 
to the forward Euler method and by (3.1) 7 we have (aa -1 )(ti) = oi(t\) = 7 so any oT 1 with 

= —1 will do. For example, the forward Euler method with a step of size — At is also an inverse 
(up to order 1 ). 

This example demonstrates that inverse methods up to order p are not unique and inverse methods 
of explicit methods need not be implicit. 



3.3 Effective order 

Effective order is achieved by using a starting method S followed by a main method M and then a 
stopping method S" 1 . We denote by a and f3 the elements of G associated with the methods M and S, 
respectively. The successive use of these three methods results in a method P = S~ 1 MS, of which the 
corresponding element of G is /3a/3 _1 . We want P to have order q, whereas M might have lower classical 
order p < q. In terms of functions in group G this leads to the following definition of the effective order 
of the Runge-Kutta method M. 

Definition 3.4. [4, Section 389] Suppose M is a Runge-Kutta method with corresponding a £ G. Then 
the method M is of effective order q if there exists a method S ( with corresponding (3 £ G) such that 

{fiafi^ 1 )^) = E(t), for every tree with r(t) < q, (3.2) 

where fj -1 is an inverse of j3 up to order q. Recall that E represents one exact step of the solution. 

The practical benefit of methods of effective order results from the observation that only method M 
need be used repeatedly, since 

P n = (S^MS) 71 = (S^MS) ■ ■ ■ (S~ 1 MS)(S- 1 MS) ~ S-^^S. 

v v ' 9 

n-times 

The starting method is applied at the beginning without advancing the solution. Instead, it introduces a 
perturbation to the solution. The main method M is then used for n time steps and finally the stopping 
method is used to correct the solution. In Section 5.2, we propose alternative starting and stopping 
procedures which allow the overall procedure to be SSP. 
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Table 3.3: Effective order five conditions of the main and starting methods M and S. We assume that 
0i =0. 

3.3.1 Effective order conditions 

For the main method M to have effective order q its coefficients must satisfy a set of algebraic conditions 
corresponding to the rooted trees of order up to q. That is, the Runge-Kutta method M corresponding 
to the function a must satisfy effective order conditions relative to the order conditions of the method S 
corresponding to the function 0. We rewrite (3.2) as (0a) (t) = (E(3)(t), for all trees with r(t) < q and 
using the product operation (3.1) we can find expressions for each tree t with r(t) < q. Each expression 
can be simplified by substituting values of ai found from previous conditions [4]. For trees up to order 
five these are tabulated in Table 3.3 (and also in [4, Sec 389]). In general, the effective order conditions 
allow more degrees of freedom on methods than the classical order conditions. Note that these conditions 
match the classical order conditions up to second order. Note also that for the tall trees ti,t4,tg,tx?, . . . 
the effective order conditions of the main method match the classical order conditions and that these are 
precisely the order conditions for linear problems [4]. 

3.4 Constructing effective order methods 

The approach we adopt is to consider the 0i as free parameters when determining the a% . The relationship 
in Table 3.3 between the ai and 0j is mostly linear (although there are a few 0| terms). It is thus 
straightforward to isolate the equations for on, determining the 0j as linear combination of the ai and 
separate the effective order conditions into conditions on the main method M and starting method S. 
This provides maximal degrees of freedom and minimizes the number of constraints when constructing 
the method M. Then when all constraints on a are found, the relative order conditions on can be 
obtained. 

The resulting effective order conditions for the main method M are given in Table 3.4 (up to effective 
order five). The order conditions for the starting method S are also given. We can also find the order 
conditions of S 1-1 in terms of the 0j (see [4, Table 386(111)]). 

Tables 3.3 and 3.4 both assume that 0i = (i.e., the starting and stopping methods perturb the 
solution but do not advance the solution in time) . This assumption is without loss of generality following 
[4, Lemma 389 A], the proof of which shows that if a method M has effective order p with particular 
starting and stopping methods (for which 0i ^ 0), then M is also effective order p with another pair of 
starting and stopping methods which do have 0i = 0. 
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Table 3.4: Effective order q, classical order p conditions on a and (3 for the main and starting methods, 
M and S respectively. 
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4 Explicit SSP Runge— Kutta methods have effective order at 
most four 

The classical order of any explicit SSP Runge-Kutta method cannot be greater than four [23]. It turns 
out that the effective order of any explicit SSP Runge-Kutta method also cannot be greater than four, 
although the proof of this result is more involved. 

Theorem 4.1. Let M denote an explicit Runge-Kutta method with C > 0. Then M has effective order 
of at most four. 

In this section we prove an even stronger result, which is stated in the following lemma. 

Lemma 4.2. Any explicit Runge-Kutta method with positive weights b > has effective order at most 
four. 

Theorem 4.1 follows immediately from Lemma 4.2 and the following well-known result. 

Lemma 4.3. (see [20, Theorem 4-2], [23, Lemma 4-2]) Any irreducible Runge-Kutta method with positive 
SSP coefficient C > must have positive weights b > 0. 

Remark 4.4. Using Lemma 4.2 and [6, Theorem 4-l]> w £ may also conclude that any explicit Runge- 
Kutta method with positive radius of circle contractivity has effective order at most four. 

The proof of Lemma 4.2 makes use of the following lemma. 

Lemma 4.5. Let b, v € M. n be given such that 

bi > for all i, (4.1a) 

n 

= 1, (4.1b) 

i=l 

/ n \ 2 

'~ ' (4-lc) 



n / n \ 2 

E hiV * = ( E biVi ) 



Then all Vi are equal but at most one; in other words, there exists fi G M and an integer k such that 
v i = A* f or a M i k- 

Proof. First observe that in the case that Vi = for all i, the stated result holds. Otherwise, let k be 
an integer between 1 and n such that Vk ^ 0. Then by collecting terms in powers of v k , (4.1c) can be 
written as 

6 fe (l - b k )v 2 k - 2b k v k biVi + J2 b i v i - ( E biV * ) = °- 

i=jtk i=jtk M^fc ' 

This is a quadratic equation in v k whose roots are real if and only if 

±b\ {V hv\ - 46 fc (l - 6*) [J2 b * v * ~ (E ) ^ °- 
Expanding and canceling terms yields 



(i-fo fe )E 6 ^ 2 - (E 6 ' v 

i^k ^i^k 



By (4.1b), 1 - b k = Y,jjtk b 3 , so we have 



E b i E biV i ~ E b o v i E biVl - °- 

j^ k i^k j^k i^k 
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Noting that the terms corresponding to i = j in the two double sums cancel and this gives 



j^k i=£k,j 

Adding the left hand side to itself, but with i,j reversed, yields 

j^k i^k,j i=£k j^k,i 

This simplifies to 

]T Y bMvi-Vj) 2 <0. 

j^k i^k,j 

Together with (4.1a), this implies that = vj for all i, j ^ k. □ 

Proof of Lemma 4-2. Any method of effective order five must have classical order at least two (see [4] 
or Table 3.4). Thus it is sufficient to show that any method with all positive weights cannot satisfy the 
conditions of effective order five and classical order two. 

Let (A, b, c) denote the coefficients of an explicit Runge-Kutta method with effective order at least 
five, classical order at least two, and positive weights b > 0. The effective order five and classical order 
two conditions (see Table 3.4) include the following: 



b T e = 1, (4.2a) 
1 

6' 



b T Ac=-, (4.2b) 



\b T c 2 - 1 = 2 , (4.2c) 
Vc 4 - b T C 2 Ac + b T (Ac) 2 = 2 2 , (4.2d) 



where the powers on vectors are understood component-wise. Define 

1 



v = -c 2 - Ac 
2 



and 

Then substituting (4.2b) in (4.2c) gives 
Also, (4. 2d) can be expressed as 



w = v 2 - (3 2 v. (4.3) 
02 = b T v. (4.4) 



/3 2 = b T v 2 . (4.5) 
Multiplying (4.4) by 02 and subtracting from (4.5) gives 

b T w = 0. (4.6) 

We divide the analysis into three cases. 



Case 1: 02 = 0. First consider the case that 02 = 0. Then b T v 2 — 0, but v / because explicit 
methods cannot have stage order two [23]. This implies that bj < for some j, which is a contradiction. 
So far we have proven the result for classical order p > 3 and the proof is similar to the result mentioned 
in [23]. The remainder of our proof deals with classical order two, where 02 ^ 0. 
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Case 2: w = 0, fa ^ 0. By the definition of w in (4.3), we have vf — favi = for all i G {1, . . . , s}, so 
for each i either Vi — or Wj = /?2- Let the set I — {i : Vi = fa}. Then (4.4) implies 

fa = = hfa = foJ2 b *> 

i=l iei iei 

which implies J2iei bi = 1> but this contradicts (4.2a) (note that v\ = because the first row of matrix 
A is identically zero). 

Case 3: w ^ 0, fa ^ 0. Since b > 0, (4.6) implies that we can choose i,j € {2, . . . , s} such that 
ifj < < Wj. By (4.3) Vi ^ 0, Uj ^ 0, and 7^ Vj. Furthermore, v\ = for any explicit method. 
Application of Lemma 4.5 reveals that all Vk must be equal except for one, which is a contradiction. 

□ 

5 Optimal effective order explicit SSP Runge— Kutta schemes 

In this section, we use the SSP theory and Butcher's theory of effective order (Sections 2 and 3) to find 
optimal explicit SSP Runge-Kutta schemes with prescribed effective order and classical order. According 
to Theorem 4.1, there are no explicit SSPRK methods of effective order five, and therefore we need only 
consider methods with effective order up to four. 

Recall from Section 3 that the methods of effective order involve a main method M as well as starting 
and stopping methods S and S^ 1 . In Section 5.2 we introduce a novel approach to construction of 
starting and stopping methods in order to allow them to be SSP. 

We denote by ESSPRK(s, q,p) an s-stage explicit SSP Runge-Kutta method of effective order q and 
classical order p. Also we write SSPRK(s, q) for an s-stage explicit SSP Runge-Kutta method of order q. 

5.1 The main method 

Our search for methods of effective order is carried out in two steps, first searching for optimal main 
methods M and then for possible corresponding methods S and S^ 1 . For a given number of stages, 
effective order, and classical order, our aim is thus to find an optimal main method, meaning one with 
the largest possible SSP coefficient C. 

To find a method ESSPRK(s, q,p) with Butcher tableau (A, b, c), we consider the optimization prob- 
lem (2.5) with &(K) representing the conditions for effective order q and classical order p (as per Ta- 
ble 3.4). The methods are found through numerical search, using Matlab's optimization toolbox. Specif- 
ically, we use fmincon with a sequential quadratic programming approach [15, 18]. This process does not 
guarantee a global minimizer, so many searches from random initial guesses are performed to help ensure 
the method with the largest possible SSP coefficient is found. 

5.1.1 Optimal SSP coefficients 

Useful bounds on the optimal SSP coefficient can be obtained by considering an important relaxation. 
In the relaxed problem, the method is required to be accurate and strong stability preserving only for 
linear, constant-coefficient initial value problems. This leads to a reduced set of order conditions and a 
relaxed absolute monotonicity condition [15, 16, 19]. We denote the maximal SSP coefficient for linear 
problems (maximized over all methods with order q and s stages) by C^ n q . 

Let C Sjq denote the maximal SSP coefficient (relevant to non-linear problems) over all methods of s 
stages with order q. Let C Sj9jP denote the object of our study, i.e. the maximal SSP coefficient (relevant 
to non-linear problems) over all methods of s stages with effective order q and classical order p. Since 
the ESSPRK(s, q,p) methods form a super class of the SSPRK(s, q) methods, we have 

Cs,q < Cs.q.p < ^s'g- (5-1) 
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9, 


123 4 5 6 7 8 9 10 11 


q = 3 


p = 2 


0.33 0.50 0.53 0.59 0.61 0.64 0.67 0.68 0.69 


q = A 


p = 2 


0.22 0.39 0.44 0.50 0.54 0.57 0.60 0.62 


q = A 


p = 3 


- - - 0.19 0.37 0.43 0.50 0.54 0.57 0.60 0.62 



Table 5.1: Effective SSP coefficients C c s — C/s of best known explicit effective order ESSPRK(s, q,p) 
methods. Entries in bold achieve the bound Cj^ given by the linear SSP coefficient and are therefore 
optimal. If no positive C can be found, we use "— " to indicate non-existence. The optimal fourth-order 
linear SSP coefficients are Cf% = 0.25 and = 0.40. 



The effective SSP coefficients for methods with up to eleven stages are shown in Table 5.1. Recall 
from Section 4 that q = 5 implies a zero SSP coefficient and from Section 3 that for q = 1, 2, the class of 
explicit Runge-Kutta methods of effective order is the simply the class of explicit Runge-Kutta methods. 
Therefore we consider only methods of effective order q — 3 and q = 4. Exact optimal values of C l " q 
are known for many classes of methods; for example see [15,16,19]. Those results and (5.1) allow us to 
determine the optimal value of Cs, g ,p a priori for the cases q = 3 (for any s) and for q — 4, s = 10, since 
in those cases we have C S:q = C x ™ q . 

5.1.2 Effective order three methods 

Since C S:q — for q = 3, the optimal effective order three methods have SSP coefficients equal to the 
corresponding optimal classical order three methods. In the cases of three and four stages, we are able 
to determine exact coefficients of families of optimal effective order methods. 

Theorem 5.1. A family of optimal three-stage, effective order three SSP Runge-Kutta methods of clas- 
sical order two, with SSP coefficient C — 1, is given by 

Yx = u n , 

Y 2 = u n + AtF{Y 1 ), 

Y 3 = u n + jAtF(Y{) + 1 AtF(Y 2 ), 

u n+i = u n + 5 ^zl A tF(Yi) + \AtF{Y 2 ) + -!-AtF(Y 3 ), 

where 1/4 < 7 < 1 is a free parameter. 

Theorem 5.2. A family of optimal four- stage, effective order three SSP Runge-Kutta methods of classical 
order two, with SSP coefficient C — 2 is given by 

Y 1 = u n , 

Y 2 = u n + ^AtF(Y 1 ), 

Y 3 = u n + -AiF(YI) + ^AtF(Y 2 ) 1 

Y A = u n + jAtF(Yi) + 1 AtF{Y 2 ) + +jAtF(Y 3 ), 



,n+l 



V; ' ^-AtF(Yx) + ~AtF(Y 2 ) + ~AtF(Y 3 ) + -^-AtF{Y 4 ), 



12 7 v 7 6 ' ' 6 v ' 12 7 
where 1/6 < 7 < 1/2 is a free parameter. 

Proof. In either theorem, feasibility can be verified by direct calculation of the conditions in problem (2.5). 



Optimality follows because C = Cl, 1 " □ 



11 



Theorem 5.1 gives a family of three-stage methods. The particular value of 7 = 1/4 corresponds 
to the classical Shu-Osher SSPRK(3,3) method [10]. Similarly, in Theorem 5.2 the particular value of 
7 = 1/6 corresponds to the usual SSPRK(4, 3) method. It seems possible that for each number of stages, 
the ESSPRK(s, 3, 2) methods may form a family in which an optimal SSPRK(s, 3) method is a particular 
member. 

5.1.3 Effective order four methods 

For effective order four, the ESSPRK(s, 4,p) methods can have classical order p = 2 or 3. In either case, 
for stages 7 < s < 11 the methods found are optimal because the SSP coefficient attains the upper bound 
of Cj™. For fewer stages, the ESSPRK methods still have SSP coefficients up to 30% larger than that of 
explicit SSPRK(s,q) methods. In the particular case of four-stage methods we have the following: 

Remark 5.3. In contrast with the non-existence of an SSPRK(4, 4) method [10,23], we are able to find 
ESSPRK(4, 4, 2) and ESSPRK(4, 4-7 3) methods. The coefficients of these methods are found in Tables 5.3 
and 5.4- 

5.2 Starting and stopping methods 

Having constructed an ESSPRK(s, q,p) scheme that can be used as the main method M, we want to find 
perturbation methods S and 5" _1 such that the Runge-Kutta scheme S~ 1 MS attains classical order q, 
equal to the effective order of method M. We also want the resulting overall process to be SSP. However 
at least one of the S and S 1-1 methods is not SSP: if j3i = then J^i &i — implies the presence of at least 
one negative weight and thus neither scheme can be SSP. Even if we consider methods with fi\ 7^ 0, one 
of S or S^ 1 must step backwards and thus that method cannot be SSP (unless we consider the downwind 
operator [9,17,24]). 

In order to overcome this problem and achieve "bona fide" effective order SSPRK methods we need to 
choose different starting and stopping methods. We consider methods R and T which each take a positive 

step such that R ~ MS and T ~ 5 I_1 M. That is, the order conditions of R and T must match those of 

9 1 

MS and 5 I_1 M, respectively, up to order q. This gives a new TM n ~ 2 R scheme which is equivalent up 
to order q to the S~ 1 M n S scheme and attains classical order q. The starting and stopping procedures 
now each take a positive step forward in time. 

To derive order conditions for the R and T methods, consider their corresponding functions in group 
G to be p and r respectively. Then the equivalence is expressed as 

p(t) = (f3a)(t) and r(t) = (a/3 -1 )(i), for all trees t with r(t) < q. (5.2) 

Rewriting the second condition in (5.2) as (r/3)(t) — a(t), the order conditions for the starting and 
stopping methods can be determined and are given in Table 5.2. These conditions could be constructed 
more generally but here we have assumed f3\ = (see Section 3.4); this will be sufficient for constructing 
SSP starting and stopping conditions. 

5.2.1 Optimizing the starting and stopping methods 

It turns out that the order conditions from (5.2) do not contradict the SSP requirements. We can thus 
find methods R and T using the optimization procedure described in Section 2.1 with the order conditions 
given by Table 5.2 for ®(K) in (2.5). 

The values of a, are determined by the main method M. Also note that for effective order q, the 
algebraic expressions on (3 up to order q — 1 are already found by the optimization procedure of the main 
method (see Table 3.4). However, the values of the order q elementary weights on (3 are not known; these 
are /J3 and /?4 for effective order three and /3s, (3q, (3? and j3g for effective order four. From Table 5.2, we 
see that both the R and T methods depend on these parameters. Our approach is to optimize for both 
methods at once: we solve a modified version of the optimization problem (2.5) where we simultaneously 
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P (t) = (M(t) 



(rP)(t)=a(t) 



Pi 


= Oil 




Tl 


= ai 




P2 


= a 2 + @2 




T2 


= a 2 




P3 


= a 3 + p 3 




7"3 


= "3 


- 2ai/3 2 - Pa 


Pi 


= «4 + a i/32 - 




T 4 


= Q?4 


~ aiP 2 - Pi 


P5 


= <*5 + Pb 




75 


= a 5 


- 3alp 2 - 3aiP 3 - p 5 


P6 


— a 6 + a 2 p 2 - 


h/9 6 


T6 


= "6 


- (a\ +a 2 ~ P 2 )P 2 - a%P 3 - aiPi - P 6 


P7 


— a 7 + ai/3 3 - 


I Pi 


7"7 


= a 7 


- 2a 1 p 4 ~ a\p 2 - p 7 


P8 


= a 8 + aiAi - 


h a 2/ 9 2 + As 


78 


= a 8 


- axP 4 - a 2 p 2 +Pi-P 8 



Table 5.2: Order conditions on p and r up to effective order four for starting and stopping methods R 
and T, respectively. The upper block represents the effective order three conditions. Here we assume 
Pi = 0. 



o 

0.730429885783319 
0.644964638145795 
1.000000000000000 



0.730429885783319 

0.251830917810810 0.393133720334985 
0.141062771617064 0.220213358584678 0.638723869798257 



0.384422161080494 0.261154113377550 0.127250689937518 0.227173035604438 
(a) Main method M, ESSPRK(4, 4, 2) 





0.545722177514735 
0.842931687441527 
0.574760809487828 
0.980872743236632 



0.545722177514735 

0.366499989048164 0.476431698393363 

0. 135697968350722 0. 176400587890242 

0.103648417776838 0.134737771331049 



0.262662253246864 
0.200625899485633 



0.541860654643112 





0.509877496215340 
0.435774135529007 
0.933203341300203 



0.233699169638954 0.294263351266422 0.065226988215286 0.176168374199685 0.230642116679654 
(b) Starting method R 



0.509877496215340 

0.182230305923759 0.253543829605247 
0.148498121305090 0.206610981494095 0.578094238501017 



0.307865440399752 0.171863794704750 0.233603236964822 0.286667527930676 
(c) Stopping method T 

Table 5.3: ESSPRK(4,4,2): an effective order four SSPRK method with four stages and classical order 
two with its associated starting and stopping methods. 



maximize both SSP coefficients subject to the constraints given in (5.2) and conditions on p given by 
Table 3.4. The unknown elementary weights on p are used as free parameters. In practice, we maximize 
the objective function min(ri, r 2 ), where n and r 2 are the radii of absolute monotonicity of the methods 
R and T. 

We were able to construct starting and stopping schemes for each main method, with an SSP coefficient 
at least as large as that of the main method. This allows the usage of a uniform time-step At < CAtpE, 
where C is the SSP coefficient of the main method. The additional computational cost of the starting and 
stopping methods is minimal: for methods R and T associated with an s-stage main method, at most 
s + 1 and s stages, respectively, appear to be required. Tables 5.3 and 5.4 show the coefficients of the 
effective SSP schemes in the case the main method is ESSPRK(4, 4, 2) and ESSPRK(4, 4, 3), respectively. 

It is important to note that in practice, if accurate values are needed at any time other than the final 
time, the computation must invoke the stopping method to obtain them. Furthermore, changing step-size 
would require first applying the stopping method with the old step-size and then applying the starting 
method with the new step-size. 
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0.601245068769724 
0.436888719886063 
0.747760163757110 



0.601245068769724 

0.139346829159954 0.297541890726109 
0.060555450075478 0.129301708677891 0.557903005003740 



0.220532078662434 0.180572397883936 0.181420582644840 0.417474940808790 
(a) Main method M, ESSPRK(4, 4, 3) 





0.438463764036947 
0.639336395725557 
0.434353425654020 
0.843416464962307 



0.438463764036947 

0.213665532574654 0.425670863150903 

0.061345094040860 0.122213530726218 

0.039559973266996 0.078812561688700 



0.250794800886942 

0.161731525131914 0.563312404874697 





0.556337718891090 
0.428870688216872 
0.815008947642716 



0.154373542967849 0.307547588471376 0.054439037790856 0.189611674483496 0.294028156286422 
(b) Starting method R 



0.556337718891090 

0. 166867537553458 0.262003150663414 
0.104422177204659 0.163956032598547 0.546630737839510 



0.203508169408374 0.096469758967330 0.321630956102914 0.378391115521382 
(c) Stopping method T 

Table 5.4: ESSPRK(4,4,3): an effective order four SSPRK method with four stages and classical order 
three with its associated starting and stopping methods. 



6 Numerical experiments 

Having constructed strong stability preserving TM n ~ 2 R schemes in the previous section, we now nu- 
merically verify their properties. Specifically, we use a convergence study to show that the procedure 
attains order of accuracy q, the effective order of M. We also demonstrate on Burgers' equation that the 
SSP coefficient accurately measures the maximal time-step for which the methods are strong stability 
preserving. 

6.1 Convergence study 

We consider the van der Pol system [12] 

u[(t) = u 2 (t), 

u' 2 (t) = - «l(t))«2(*) - Ul(t), 



over the time interval t £ [0,50] with fj, = 2 and initial values tii(0) = 2 and 1*2(0) = 1. The reference 
solution for the convergence study is calculated by Matlab's ode45 solver with relative and absolute 
tolerances set to 10~ 13 . 

We solve the initial value problem (6.1) using SSP TM n ~ 2 R schemes. The solution is computed 
using n — 100 • 2 k time steps for k = 2, . . . , 7. The error at t — 50 with respect to time-step is shown 
in Figure 6.1 on a logarithmic scale. The convergence study is performed for TM n ~ 2 R schemes with 
various number of stages s and the results show that the schemes attain an order of accuracy equal to 
the effective order of their main method M . It is important in doing this sort of convergence study that 
the effective order can only be obtained after the stopping method is applied. Intermediate steps will 
typically only be order p accurate (the classical order of the main method). Finally, we note that the 
methods with more stages generally exhibit smaller errors (for a given step size). 
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10 10 10 10 10 10 10 10 

At At 
(a) ESSPRK(s, 3, 2) (b) ESSPRK(s, 4, 2) 

Figure 6.1: Convergence study of TM n ~ 2 R Runge-Kutta schemes when (a) M is an ESSPRK(s, 3, 2) 
method and (b) M is an ESSPRK(s, 4, 2) method. 

6.2 Burgers' equation 

The inviscid Burgers' equation consists of the scalar hyperbolic conservation law 

U t + f(U) x =0, (6.2) 

with flux function f(U) — \U 2 . We consider initial data U{0,x) = \ — |sin7ra, on a periodic domain 
x G [0,2). The solution advances to the right where it eventually exhibits a shock. We perform a 
semi-discretization using an upwind approximation to obtain the system of ODEs 

d f(ui)-f(ui-i) 
At Ax 

This spatial discretization is total-variation-diminishing (TVD) when coupled with the forward Euler 
method under the restriction [21] 

At < At FE = Az/ 1| f/^aOHoo. 

Recall that a time discretization with SSP coefficient C will give a TVD solution for At < CAtpE- 

Burgers' equation was solved using an SSP TM n ~ 2 R scheme with time-step restriction At < crAtpE, 
where a indicates the size of the time step. We integrate to roughly time tf = 1.62 with 200 points in 
space. Figure 6.2 shows that if a is chosen less than the SSP coefficient of the main method, then no 
oscillations are observed. If this stability limit is violated, then oscillations may appear. 

We measure these oscillations by computing the total variation of the numerical solution. When M 
is an ESSPRK(4, 4, 2) method, it turns out that a = 1.57 is the largest value of a for which the total 
variation is monotonically decreasing during the calculation. This is 79% larger than the value of the 
SSP coefficient C = 0.88. 

We also consider Burgers' equation with a discontinuous square wave initial condition 

TTta \ / !' 0.5 < a; < 1.5 , , 

(7(0,2:) = < n . (6.3) 

[0, otherwise. 

The solution consists of a rarefaction (i.e., an expansion fan) and a moving shock. Again we use 200 
points in space and we compute the solution until roughly time tf = 0.6. Figure 6.3 shows the result 
of solving the discontinuous problem using an SSP TM n ~ 2 R scheme, where M is an ESSPRK(5, 4, 2) 
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i = 1.6253 



TV = 0.96028 



t = 1.6213 



TV = 1.0075 




(a) a = 0.88 




0.6 

X 



(b) cr = 1.60 



Figure 6.2: Solution of Burgers' equation at the final time with continuous initial data, using a TM n 2 R 
scheme, where M is the optimal ESSPRK(4, 4, 2). The SSP coefficient is C = 0.88. Figure 6.2b shows a 
zoom in the region of space where oscillations appear. Here TV denotes the TF-norm of the solution at 
the final time: a value greater than 1 (the TV-norm of the initial condition) indicates a violation of the 
TVD condition. 



q, p~~^^ 


3 4 56789 10 11 


9 = 3 


p = 2 


1.04(4%) 2.00(0%) 2.65(0%) 3.52(0%) 4.29(0%) 5.11(0%) 6.00(0%) 6.79(0%) 7.63(0%) 


q = A 


p = 2 


1.07(22%) 1.98(2%) 2.69(2%) 3.56(1%) 4.33(1%) 5.16(1%) 6.05(1%) 6.84(1%) 


9 = 4 


p = 3 


1.05(35%) 1.89(3%) 2.63(2%) 3.53(1%) 4.31(1%) 5.16(1%) 6.04(1%) 6.85(1%) 



Table 6.1: Maximum observed coefficients exhibiting the TVD property on the Burgers' equation ex- 
ample with discontinuous data (6.3). The numbers in parenthesis indicate the increase relative to the 
corresponding SSP coefficients. 



method with SSP coefficient C = 1.95. In this case, a = 1.98 appears to be the largest value for which 
the total variation is monotonically decreasing during the calculation. This is 2% larger than the value of 
the SSP coefficient. Figure 6.3b shows part of the solution exhibiting oscillations when a is larger than 
the SSP coefficient. For various schemes, Table 6.1 shows the maximum observed values of a for which 
the numerical solution is total variation decreasing for the entire computation. With the exception of the 
four-stage effective order four methods, we note good agreement between the SSP coefficient predicted 
by the theory and the maximum time-step for which the numerical solution is TVD. 

We also note the necessity of our modified starting and stopping methods in the RM n ~ 2 T approach: 
in this example if we use the original approach of S and S , the solution exhibits oscillations immediately 
following the application of the starting perturbation method S. 



7 Conclusions 

We use the theory of strong stability preserving time discretizations with Butcher's algebraic interpreta- 
tion of order to construct effective order SSP Runge-Kutta (ESSPRK) methods. These methods, when 
accompanied by starting and stopping methods, attain an order of accuracy higher than their (classical) 
order. We propose a new choice of starting and stopping methods to allow the overall procedure to be 
SSP. We prove that explicit Runge-Kutta methods with strictly positive weights have at most effective 
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t = 0.60554 TV = 2 ( = 0.602 TV = 2.0006 




0.5 1 1.5 2 1.3 1.4 1.5 1.6 1.7 1.8 



X X 
(a) a = 1.95 (b) a = 2.15 

Figure 6.3: Solution of Burgers' equation at the final time with discontinuous initial data, using a 
TM n ~ 2 R scheme, where M is ESSPRK(5, 4, 2) method. The SSP coefficient is C = 1.95. Figure 6.3b 
shows a zoom in the region of space where oscillations appear. Here TV denotes the TV-norm of the 
solution at the final time: a value greater than 2 indicates a violation of the TVD condition. 

order four. This extends the barrier already known in the case of classical order explicit SSPRK methods. 

ESSPRK methods of effective order three and four are constructed by numerical optimization. Most 
of the methods found are optimal because they achieve the upper bound on the SSP coefficient known 
from linear problems. Also, despite the non-existence of four-stage, order four explicit SSPRK methods, 
we find effective order four methods with four stages (of classical order two and three). We perform 
numerical tests which confirm the accuracy and SSP properties of the ESSPRK methods. 

The ideas here are applied to explicit Runge-Kutta methods, but they could also be applied to other 
classes of methods including implicit Runge-Kutta methods, general linear methods, and Rosenbrock 
methods. 
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